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Interpolation  approximations 
for  M\G\oo  arrival  processes 
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Abstract 


We  present  an  approximate  analysis  of  a  discrete-time  queue  with  correlated  ar¬ 
rival  processes  of  the  so-called  M|G|oo  type.  The  proposed  heuristic  approximations 
are  developed  around  asymptotic  results  in  the  heavy  and  light  traffic  regimes.  Inves¬ 
tigation  of  the  system  behavior  in  light  traffic  quantifies  the  differences  between  the 
gradual  M\G\oo  inputs  and  the  point  arrivals  of  a  classical  GI\GI\l  queue.  In  heavy 
traffic,  salient  features  are  effectively  captured  by  the  exponential  distribution  and  the 
Mittag-Leffler  special  function,  under  short-  and  long-range  dependence  respectively. 
By  interpolating  between  the  heavy  and  light  traffic  extremes  we  derive  approxima¬ 
tions  to  the  queue  size  distribution,  applicable  to  all  traffic  intensities.  We  examine  the 
accuracy  of  these  expressions  and  discuss  possible  extensions  of  our  results  in  several 
numerical  examples. 


1  Introduction 

The  conclusions  of  a  series  of  measurement  studies  demonstrating  that  network  traffic  ex¬ 
hibits  persistent  long  term  correlations  have  spurred  recent  activity  in  the  study  of  queueing 
systems  with  correlated  arrival  processes.  Analytical  results  reveal  that,  when  strong  de¬ 
pendencies  are  present,  diverse  queueing  patterns  may  arise,  in  contrast  to  the  familiar 
exponential  decay  encountered  in  traditional  traffic  models  with  bounded  exponential  mo¬ 
ments. 

*The  work  of  this  author  was  supported  through  NSF  Grant  NSFD  CDR-88-03012  and  the  Army  Research 
Laboratory  under  Cooperative  Agreement  No.  DAAL01-96-2-0002. 

iThe  work  of  this  author  was  supported  partially  through  NSF  Grant  NSFD  CDR-88-03012,  NASA  Grant 
NAGW277S  and  the  Army  Research  Laboratory  under  Cooperative  Agreement  No.  DAAL01-96-2-0002. 
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In  this  paper  we  deal  with  a  discrete  time  queue,  viewed  as  a  surrogate  for  a  network 
multiplexer,  driven  by  an  M | G | oc  arrival  stream.  Both  the  discrete  time  M\G\oo  process 
considered  here  and  its  continuous  time  variant  are  among  the  traffic  models  arising  from 
large  aggregations  of  on/off  sources,  that  have  attracted  a  great  deal  of  attention.  Reasons 
for  this,  such  as  flexibility  in  representing  correlation  functions  of  actual  traffic  traces,  are 
discussed  in  [9].  A  fluid  queue  fed  by  the  continuous  time  version  of  the  process  has  been 
studied  at  least  as  early  as  1974  [3].  Later,  Cox  notices  that  the  M\G\oo  busy  server  process 
with  heavy  tailed  G  is  a  second  order  asymptotically  self-similar  process  [4], 

Here,  we  are  interested  in  the  entire  steady-state  queue  length  distribution  at  a  multi¬ 
plexer  fed  by  the  M\G\oo  arrival  process.  For  arbitrary  pmf  G.  the  system  lacks  the  desired 
Markovian  structure  and  a  calculation  using  numerical  inversion  techniques,  e.g.  [13],  is  not 
possible,  since  no  ^-transform  expressions  are  available.  An  exception  is  provided  by  the  ge¬ 
ometric  case,  where  a  two  dimensional  Markov  chain  formulation  and  a  functional  equation 
for  the  ^-transform  is  given  in  [2], 

To  circumvent  the  difficulties  of  an  exact  analysis,  one  may  rely  on  information  gleaned 
from  various  asymptotic  regimes.  A  promising  approach  consists  of  deriving  approximations 
from  the  analysis  of  large  buffer  asymptotics  [7,  8,  10];  these  estimates  are  exact  in  the  limit 
as  the  buffer  level  goes  to  infinity. 

Our  objective  is  to  explore  alternative  approximations  to  the  queue  length  probabilities, 
developed  around  a  combination  of  fight  and  heavy  traffic  asymptotics.  Such  approximations 
become  exact  in  the  limit  as  the  traffic  intensity  goes  to  zero  and  one  respectively.  In  fight 
traffic  we  take  advantage  of  the  fact  that  the  M\G\oo  arrival  process  is  obviously  “Poisson 
driven”,  so  that  the  Reiman-Simon  theory  [12]  applies,  under  a  bounded  exponential  moment 
assumption.  The  resulting  light  traffic  limits  of  the  queue  with  M\G\oo  arrivals  differ  from 
those  of  a  classical  GI\GI\1  queue.  This  is  a  manifestation  of  the  fact  that  work  that  joins 
the  system  gradually,  as  is  the  case  with  M\G\oo  inputs,  generates  less  queueing  than  work 
that  arrives  instantaneously.  From  the  heavy  traffic  regime,  we  collect  the  associated  limit 
distribution  of  the  queue  size:  This  is  given  through  the  exponential  function  in  the  standard 
short-range  dependent  setup,  and  the  Mittag-Leffler  special  function  in  the  case  where  the 
M\G\oo  process  is  long-range  dependent.  The  approximation  to  the  queue  size  distribution 
is  subsequently  generated  by  interpolating  between  the  heavy  and  fight  traffic  extremes.  For 
some  common  pmfs  G  the  approximant  assumes  a  simple  final  form.  More  interestingly,  it 
has  the  potential  of  capturing  accurately  the  queue  size  distribution  at  small  buffer  sizes, 
for  which  approximations  based  on  large  buffer  asymptotics  are  usually  ill  fitted.  On  the 
other  hand,  when  G  has  finite  exponential  moment,  we  do  not  expect  the  heavy-light  traffic 
interpolation  to  be  accurate  for  buffer  sizes  much  larger  than  the  maximum  burst  length:  It 
simply  does  not  possess  the  correct  decay  rate  -  it  does  so  only  as  the  traffic  intensity  tends 
to  one,  i.e. ,  in  the  heavy  traffic  limit.  Surprisingly,  this  drawback  is  often  absent  under  long- 
range  dependence,  since  there  are  cases  where  the  queue  size  distribution  has  hyperbolic 
asymptotics  with  the  same  exponent  for  all  traffic  intensities!  Then  an  approximation  is 
more  valuable,  especially  when  considering  that,  in  the  presence  of  heavy  tails,  alternative 
estimates  by  means  of  simulation  take  an  unreasonably  long  time  to  obtain.  Yet,  this  is 
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somewhat  compromised  by  the  unaivalability  of  rigorously  established  light  traffic  limits, 
under  long-range  dependence.  In  Section  6  we  rely  on  a  postulated  relationship,  but  the 
problem  is  still  unresolved. 

The  paper  is  organized  as  follows:  The  description  of  the  system  model  is  given  in  Section 
2.  Section  3  contains  the  main  conclusions  of  the  light  and  heavy  traffic  analyses.  These  are 
the  ingredients  for  the  approximation,  which  is  presented  in  Section  4  and  discussed  through 
numerical  examples  in  Section  5.  Further  extensions  of  the  results  are  suggested  in  Section 
6. 


2  The  System  Model 


We  introduce  the  queueing  model  of  interest,  together  with  the  required  notation.  We  start 
by  presenting  the  M\G\oo  arrival  processes  and  several  of  their  properties;  additional  facts 
can  be  found  in  [4,  10]. 

Consider  a  population  of  infinitely  many  information  sources,  operating  in  discrete-time. 
Sources  can  be  in  one  of  two  states,  active  or  idle.  During  time  slot  [n,  n  +  1),  n  =  0, 1, . . ., 
Pn+ 1  new  sources  become  active.  Source  j,  j  =  1, . . .  ,(3n+i  begins  generating  information  by 
the  start  of  slot  [n+ 1,  n+2),  its  activity  period  has  duration  crn+i,j  (in  number  of  slots).  While 
active,  each  source  emits  information  at  a  constant  rate  of  one  information  unit  (packet)  per 
time  slot.  After  its  activity  period  expires,  each  source  switches  off  permanently,  never  to 
generate  packets  again.  Let  bn  denote  the  number  of  active  sources,  or  equivalently,  the 
number  of  packets  generated  by  the  active  sources  at  the  beginning  of  time  slot  [n,  n  +  1). 
If  initially  (i.e.,  at  time  n  =  0)  there  were  already  b  active  sources,  we  denote  by  cr0j  the 
residual  activity  duration  (in  time  slots)  for  the  jth  active  source,  j  =  1 , ,b. 

Throughout,  the  FV-v alued  rvs  b,  {/3n+i,  n  =  0, 1, . . .},  {<Jn,j,  n  =  1,2,...;  j  —  1,2,...} 
and  {(J0j,  j  =  1,2,.. .}  satisfy  the  following  assumptions:  (i)  These  rvs  are  mutually  inde¬ 
pendent;  (ii)  The  rvs  {/?n+i,  n  —  0, 1, . . .}  are  i.i.d.  Poisson  rvs  with  parameter  A  >  0;  (iii) 
The  rvs  {crn,j,  n  =  1,  -  -  - ;  j  =  1,2,...}  are  i.i.d.  with  common  pmf  G  on  {1,  2, . . .}.  Let 
a  be  a  generic  IN  valued  rv  distributed  according  to  the  pmf  G,  assume  throughout  that 
E  [cr]  <  cx);  (iv)  The  rvs  {cr0j,  j  =  1,2,...}  are  i.i.d.  FV-valued  rvs  distributed  accord¬ 
ing  to  the  equilibrium  pmf  Ge  associated  with  G,  i.e.,  if  ae  denotes  a  generic  FV-valued  rv 
distributed  according  to  the  pmf  Ge,  then 

P  \ae  =  n]  =  P  ~  U ^ ,  n  =  l,2,...  (1) 

hi  [(j 


In  summary,  the  process  {bn,  n  =  0, 1, . . .}  results  from  discrete-time  Poisson(A)  arrivals 
of  information  sessions,  where  the  session  duration  is  distributed  according  to  the  pmf  G 
and  the  packet  generation  rate  of  an  ongoing  session  is  one  packet  per  time  slot.  Under  the 
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enforced  assumptions  {bn,  n  =  0, 1, . . .}  can  be  identified  as  the  busy  server  process  of  a 
discrete-time  M\G |oo  queue;  for  this  reason  the  packet  arrival  process  { bn ,  n  =  0, 1, . . .}  is 
referred  to  as  the  M\G\oo  arrival  process.  The  following  proposition  shows  that  { bn .  n  = 
0, 1, . . .}  is  a  correlated  process,  with  time  dependencies  controlled  by  the  tail  of  a  [10]. 


Proposition  2.1  If  b  is  taken  to  be  a  Poisson  rv  with  parameter  AE  [cr],  then  the  process 
{bn,  n  =  0, 1, . . .}  is  a  (strictly)  stationary  ergodic  process  with  the  properties: 

(a)  For  each  n  =  0, 1, . . .,  the  rv  bn  is  a  Poisson  rv  with  parameter  AE  [cr]  ; 

(b)  Its  covariance  function  is  given  by 


co  v(bn+j,bn) 


AE 


J']  + 


AE  [a]  P  [ae  >  j] ,  n,j  =  0,1,... 


(c)  Its  index  of  dispersion  of  counts  (IDC)  is  given  by 


oo  OO  \ 

IDC  =  Y  cov ( bn+j ,  bn)  =  AE  [a]  Y  p  [ae  >  j)  =  -E  [cr(cr  +  1)]  , 

3=0  3=0  1 


and  the  process  is  short-range  dependent  (i.e.,  IDC  finite)  if  and  only  if  E  [cr2]  is  finite. 


We  now  feed  this  M G  |  oc  arrival  stream  { bn ,  n  —  0, 1, . . .}  into  a  discrete-time  single 
server  queue  with  infinite  buffer  capacity.  Such  a  queueing  system  routinely  serves  as  a 
model  for  a  network  multiplexer:  If  qn  denotes  the  number  of  packets  remaining  in  the 
multiplexer  buffer  by  the  end  of  slot  [n  —  1,  n),  and  the  multiplexer  output  link  can  transmit 
c  packets/slot,  then  the  buffer  content  sequence  { qn ,  n  —  0, 1, . . .}  evolves  according  to  the 
Lindley  recursion 

q0  =  0;  qn+i  =  [q„  +  bn+i  -  c]  +  ,  n  =  0, 1, . . .  (2) 

From  Part  (a)  of  Proposition  2.1  the  average  input  rate  to  the  multiplexer  is  E  [bri]  = 
AE  [cr],  and  the  system  is  stable  if  the  traffic  intensity  p  =  AE  [cr]  /c  satisfies  p  <  1.  In  that 
case  qn  =^n  q,  where  the  IR-valued  rv  q  is  the  stationary  queue  size  in  the  multiplexer 
buffer.  We  are  interested  in  evaluating 

P(b,  p)  =  Pp  [q  >  b] ,  b>  0  and  Qm(p )  =  Ep  [q'n] ,  0  <  p  <  1,  m  =  1,2, 

i.e.,  the  probability  that  the  stationary  queue  size  exceeds  b,  and  the  queue  size  first  and  sec¬ 
ond  moments,  when  the  traffic  intensity  is  p.  To  that  end  we  develop  simple  approximations 
that  all  flow  from  asymptotic  results  under  heavy  and  light  traffic  conditions. 


3  Heavy  and  Light  Traffic 


The  interpolation  approximation  we  have  in  mind  hinges  on  the  availability  of  explicit  ex¬ 
pressions  for  limits  of  system  quantities  as  p  — »  1  (heavy  traffic  limits),  and  derivatives  with 
respect  to  p  as  p  — »  0  (light  traffic  derivatives).  It  thus  requires  examination  of  the  behavior 
of  the  queue  with  M\G\oo  arrivals  under  each  one  of  these  two  asymptotic  regimes. 
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Light  Traffic  We  start  with  the  light  traffic  regime.  The  right-hand  derivatives  at  p  =  0 
of  the  various  metrics  of  interest  are  evaluated  using  the  Reiman-Simon  technique  [12].  For 
the  system  to  be  in  the  domain  of  applicability  of  the  Reiman-Simon  results,  an  assumption 
on  finiteness  of  the  exponential  moment  of  o  is  needed. 


Assumption  (A) 


There  exists  9*  >  0  such  that  E 


<  oo  for  9  <  9*. 


The  detailed  fight  traffic  analysis  of  the  queue  with  M\G\oo  arrivals  is  provided  in  [15].  We 
summarize  the  conclusions  in  the  following. 


Proposition  3.1  Consider  the  Lindley  recursion  (2)  with  integer  release  rate  c  =  1,2,.. 
and  let  b  =  0,1,  ... .  Under  Assumption  (A)  it  holds  that 

(a)  For  each  n  =  0, 1, . . . ,  c 

J^P(&,0+)  =  0  and  Qm(  0+)  =  0,  m  =  l,2, ....  (3) 

(b)  In  addition,  for  c  =  1 , 


,2  <92 


dp' 


and 


=  E 


2  <3i(o+) 


ip  °2(0+) 


a  -  b}+2]  P  [a  >  b]  +  2  E  f  [cr  —  6] 


1  2 


— 3  E 

Ek!l 

EM 


a  —  b]+  P  [cr  >  6]  +  P  [cr  >  bf 


1  |  E[cr2]2' 

2  l1+  EH2, 


(4) 

(5) 

(6) 


Proposition  (3.1)  delineates  a  light  traffic  behavior  for  the  queue  with  M\G\oo  arrivals 
that  is  certainly  different  from  the  one  of  a  classical  GI\GI\1  queue.  As  seen  from  (3),  when 
c  =  1,  in  which  case  the  multiplexer  can  serve  no  more  than  one  source  per  time  slot,  the  first 
derivative  of  the  tail  probability  is  zero.  Hence,  in  a  Taylor  expansion  of  Pp  [q  >  b]  around 
p  =  0  the  linear  term  in  p  would  offer  no  contribution.  In  contrast,  the  stationary  workload 
W  in  a  single  server  M\G\1  queue  is  known  to  satisfy  Pp  [W  >  x]  ~  p(l  —  Ge(x ))  (p  — *  0), 
that  is,  in  the  classical  queueing  setup  the  corresponding  expansion  starts  with  a  non-zero  p 
term.  For  M\G\oo  arrivals  it  is  the  second  derivative  (4)  which  is  the  most  informative.  This 
highlights  the  role  of  the  activity  duration  rv  a,  through  both  its  distribution  and  its  first 
two  moments.  Notice  that  even  if  Assumption  (A)  were  to  be  relaxed,  (4)  shows  that  for 
Pp  [ q  >  b]  to  decay  like  p2  for  small  p  it  is  necessary  that  E  [a2]  be  finite.  If  E  [a2]  =  oo,  as 
is  the  case  for  long-range  dependent  M\G\oo  arrivals,  expression  (4)  yields  infinity  and  p2  is 
no  longer  the  correct  order  of  decay.  A  different,  perhaps  smaller  exponent  should  be  sought 
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in  the  long-range  dependent  case  (see  (21)).  Finally,  relations  (3)  reflect  (though  in  a  rough 
manner)  the  statistical  multiplexing  gain:  Since  the  first  non-zero  contribution  to  a  the  tail 
probability  is  no  lower  than  pc+1,  (3)  implies  that  increasing  the  multiplexer  release  rate  c 
while  maintaining  the  same  traffic  intensity  p  would  result  in  a  decreasing  tail  probability 
Pp  [q  >  b],  as  could  be  expected. 


Heavy  traffic  Considered  next  is  the  behavior  of  the  queue  with  M\G\oo  arrivals  in  heavy 
traffic,  that  is,  as  the  arrival  rate  AE  [cr]  tends  to  the  multiplexer  release  rate  c  from  below. 
Clearly,  as  the  traffic  intensity  p  converges  to  one  the  system  becomes  unstable  and  the 
queue  length  grows  unbounded.  It  is  thus  necessary  to  seek  a  suitable  normalizer  for  the 
queue  length  process,  so  that  its  normalized  version  has  a  non  trivial  heavy  traffic  limit. 
This  problem  is  typically  addressed  in  a  setup  where  the  system  of  interest  is  embedded  into 
a  family  of  queueing  systems,  parametrized  by  an  integer,  say  l  =  1,2,...,  ensuring  that,  as 
l  t  oo,  the  appropriate  trend  to  instability  is  established.  Such  an  approach  was  pursued 
in  [16],  providing  a  complete  characterization  of  the  arising  heavy  traffic  limits.  We  tacitly 
assume  here  that  the  heavy  traffic  limit  of  the  stationary  distribution  coincides  with  the 
stationary  distribution  of  the  heavy  traffic  limit,  and  gather  the  required  results  from  [16]  in 
a  convenient  form: 


Proposition  3.2  The  heavy  traffic  limits  of  the  stationary  queue  length  distribution  associ¬ 
ated  with  (2)  can  be  classified  as  follows: 

(a)  If  E  [a2]  <  oo  then 


limPp  [(1  —  p)  q  >  x]  =  exp 

p-* 1 


2E  [a 

'  E  [cr2] 


X 


x  >  0. 


(b)  If  P  [a  >  n\  =  n  ",  n  =  1,2,...,  with  1  <  a  <  2,  then 


liniP, 

p->i  * 


(1  ~P) 


!/(«-!) 


q  >  x 


—  Ea- 1 


*>o, 


where 


Eu(x)  = 


x 


Wo  F(^n  +  1) 
is  the  Mittag-Leffler  special  function  [5]. 


V  r(2-a) 

,  v  >  0,  x  G  R, 


(7) 

(8) 

(9) 


Part  (a)  of  Proposition  3.2  addresses  the  classical  short-range  dependent  case,  for  which 
the  heavy  traffic  normalizer  is  (1 —p)  and  the  limiting  heavy  traffic  distribution  is  exponential. 
Part  (b)  deals  with  a  long-range  dependent  M\G\oo  arrival  process.  Under  long-range 
dependence,  the  heavy  traffic  queue  length  distribution  is  expressed  through  a  Mittag-Leffler 
function  with  hyperbolic  decay,  while  the  power-law  behavior  of  the  heavy  traffic  normalizer 
is  (1  —  p)1^"-1).  We  mention  that  this  result  can  be  stated  in  a  more  general  manner  to 
cover  the  situation  where  the  tail  of  a  is  regularly  varying  of  order  a,  1  <  a  <  2. 
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The  results  under  the  light  and  heavy  traffic  regimes  are  subsequently  combined  into 
approximations  for  all  values  of  the  traffic  intensity. 


4  Interpolation  Approximations 


Whenever  Assumption  (A)  is  satisfied,  Pp  [q  >  b]  is  infinitely  differentiable  with  respect  to 
p  at  p  =  0,  hence  it  can  be  approximated  by  bringing  together  heavy  traffic  limits  and  light 
traffic  derivatives  into  a  Taylor  series-like  expansion.  To  this  end  we  follow  the  approach 
proposed  in  [6].  In  passing,  we  also  discuss  approximations  for  Qm(p),  m  —  1,2.  The  details 
are  as  follows: 


The  interpolation  method  Consider  the  normalized  queue  length  rv  (1  —  p)  q  and  define 
F(x ,  p)  =  Pp  [(1  —  p)  q  >  x] ,  0  <  p  <  1,  x  >  0  (10) 

and 

F(x,  1)  =  limPp  [(1  —  p)  q  >  x]  .  (11) 

Assume  that  partial  derivatives  of  F(x,p)  with  respect  to  p.  up  to  order  n.  at  p  =  0+, 
are  available.  Construct  Fn(x,p),  the  nth  order  interpolation  approximation  to  F(x,p),  by 
means  of  the  polynomial 

Ft  i,  /  77/1  C\%  \ 

F»(x,  P)=  F(X,  0+)  +  (f(z,  1)  -  S  F(x,  0+))  p“+1.  (12) 


Observe  that 


Fn(x,  1)  =  F(x,  1)  and  Fn(x,  0+)  =  F(x,  0+),  i  =  0,l,...,n, 

that  is,  Fn(x,p)  is  precisely  that  unique  n  +  1  degree  polynomial  in  p  which  matches  the 
n  +  1  partial  derivatives  of  F(x,  p)  at  p  —  0+  and  its  heavy  traffic  limit.  Now,  by  reversing 
the  (1  —  p)  normalization  in  Fn(x,p)  we  generate  the  nth  order  interpolation  approximation 
to  Pp  [q  >  b]  as 

Pp[tf  >  b\  «  F„((l  ~  p)  b,p).  (13) 

Note  that,  in  principle,  this  may  lie  outside  [0,1],  in  which  case  it  is  obviously  a  poor 
approximation.  To  calculate  the  quantities  associated  with  (13)  it  remains  to  express  the 
partial  derivatives  appearing  in  (12)  in  terms  of  the  light  traffic  derivatives  of  Pp  [q  >  b\ .  We 
have 

^  F(x,0+)  =  -^  P(x,0+)  +  x  ^P(x,0+)  (14) 
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and 


dp 2 


F(x,  0+)  = 


d 2  d2 

P(x,  0+)  +  2x  P(x,  0+) 


dp 


dpdx 


+2x  —  P(x,  0+)  +  x2  - t-^P(x,  0+). 
dx  dxz 


(15) 


In  case  additional  light  traffic  information  is  available,  repeated  application  of  the  chain  rule 
will  yield  higher  order  derivatives,  as  needed. 


Approximate  expressions  We  are  now  ready  to  write  approximate  expressions  anchored 
on  the  heavy  and  light  traffic  results  of  Section  3.  Proposition  3.2(a)  provides  the  limit  (11) 
that  should  be  inserted  in  (12).  Proposition  3.1(a)  can  be  used  to  substitute  for  the  partials 
in  (14)  and  then  in  (12).  Thus,  if  the  multiplexer  release  rate  is  c  =  1,2, . . .,  the  cth  order 
interpolation  approximation  to  Pp  [q  >  b]  is  simply 

p p  [q  >  b]  ~  Fc((  1  -  p)  b ,  p)  =  pc+l  exp  (!  ~  p)  ^  •  (16) 

More  can  be  accomplished  in  the  case  c  =  1,  since  Proposition  3.1(b)  affords  us  a  promising 
2nd  order  interpolation  approximation.  Insertion  of  (15)  in  (12)  yields 

Hb,  p)  =  ^p2(l  ~  p)  P(b,  0+)  +  p 3  exp  (-frpj  (17) 

and  the  latter  leads  to  the  2nd  order  approximant 

Pp[<7  >  b]  w  F2((l  -  p)  b,p),  c=  1,  (18) 

where  Proposition  3.1(b)  is  used  to  supply  the  second  partial  derivative  in  (17). 


Next,  we  briefly  deal  with  moment  approximations.  We  restrict  attention  to  the  case 
c  =  1  and  consider  only  the  queue  length  first  and  second  moment.  The  relevant  light  traffic 
limits  are  given  by  (3),  (5)  and  (6).  In  heavy  traffic  it  can  be  inferred  from  (7)  that 

lim(l  -  p)m  Qm(p )  =  ml  ,  m  —  1,2, ... . 


Moment  interpolations  are  then  developed  in  very  much  the  same  manner  as  distribution 
interpolations.  We  skip  the  details  of  the  derivation  and  list  the  final  expressions 


and 


Qi(p) 


E  \c2\  p2 
2E\a\  1  -  p 


c  =  1 


(19) 


Qi{p) 


'Ek 


21 2 


4(1  —  p)2 


EM 


1  + 


Eb 


2p 


EM 


+  1 


c  =  1. 


(20) 


We  stress  that  formula  (19)  is  in  fact  exact.  This  is  a  result  whose  continuous  time  analog  has 
been  established  for  a  more  general  fluid  model  in  [14,  p.  23].  On  the  contrary  (20)  cannot 
be  exact.  To  verify  this  consider  the  example  where  a  =  1  deterministic.  This  corresponds 
to  i.i.d.  Poisson  arrivals,  for  which  a  probability  generating  function  of  the  queue  length  rv 
is  available.  It  can  be  shown  [15]  that  the  exact  expression  is 

P2 

Q 2(p)  =  (p2  —  P  +  3),  C=l,  CT  =  1  a.s. 

a  formula  that  clearly  cannot  be  recovered  using  only  two  light  traffic  derivatives.  Still,  when 
a  =  1  approximation  (20)  is  within  9%  of  the  correct  value,  for  all  traffic  intensities. 


We  close  this  section  with  a  comment.  Recall  that  each  active  source  in  the  M | G | oc 
arrival  process  generates  one  information  unit  per  time  slot.  So,  c  =  1  corresponds  to  the 
case  where  the  amount  of  service  in  one  slot  is  exactly  equal  to  the  amount  of  information 
that  one  active  source  generates  in  one  slot.  When  c  =  1  a  single  active  source  suffices  to 
make  full  use  of  the  server  capacity;  in  this  system  there  is  never  any  leftover  capacity  to 
simultaneously  serve  more  than  one  sources.  On  the  contrary,  when  c  >  1,  the  server  can 
attend  to  more  than  one  sources  during  one  time  slot,  so  that  there  is  a  multiple  service 
feature  to  the  system  behavior.  An  exact  or  approximate  analysis  in  this  regime  is  clearly 
more  challenging. 


5  Numerical  Results 


To  gauge  the  accuracy  of  the  proposed  expressions  we  carry  out  simulation  experiments 
choosing  various  distributions  for  the  burst  duration  rv  a.  The  experimental  values  are 
obtained  by  regenerative  simulation  and  relative  widths  accompanying  them  correspond  to 
95%  confidence  intervals.  We  confine  ourselves  to  the  simple  situation  where  the  multiplexer 
release  rate  is  c  =  1.  While  the  list  of  examples  below  is  not  exhaustive,  it  does  serve 
to  illustrate  the  ability  of  the  heavy-light  traffic  interpolation  to  “ballpark”  the  true  tail 
probabilities,  as  well  as  its  limitations. 


Deterministic  When  the  burst  duration  is  deterministic,  a  =  D  a.s.,  D  =  1,2,..., 
approximation  (18)  reads 

p ',[«>&]  «  dAld  (3[£>  -  (1  -  p)6]+([D  -  (1  -  p)b]+  -  1) 

+  1  [D  >  (1  —  p) 6] )  +  p3 exp  (-4(1  -  p)Pj  ,  b  =  0,1,.... 

We  let  the  burst  duration  be  cr  =  3  and  obtain  simulation  estimates  for  the  steady  state 
probability  Pp  [q^  >  0].  In  this  case  an  explicit  expression  for  Pp  \q^  >  0]  is  available  [15]. 
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Tail  probability  P  [q^  >  0] 

p 

Exact 

Simulation 

Approximation 

Error  (%) 

0.1 

1.0478e-02 

1.0469e-02±0.2% 

1.0500e-02 

-0.21 

0.2 

4.1622e-02 

4.1668e-02±0.3% 

4.1778e-02 

-0.38 

0.3 

9.3042e-02 

9.3020e-02±0.2% 

9.3500e-02 

-0.49 

0.4 

1.6441e-01 

1.6442e-01±0.2% 

1.6533e-01 

-0.56 

0.5 

2.5545e-01 

2.5533e-01±0.2% 

2.5694e-01 

-0.58 

0.6 

3.6594e-01 

3.6607e-01±0.1% 

3.6800e-01 

-0.56 

0.7 

4.9573e-01 

4.9601e-01±0.1% 

4.9817e-01 

-0.49 

0.8 

6.4470e-01 

6.4488e-01±0.1% 

6.4711e-01 

-0.37 

0.9 

8.1279e-01 

8.1264e-01±0.1% 

8.1450e-01 

-0.21 

Table  1:  P  [q^  >  0]  for  deterministic  burst  duration  a  =  3. 

In  Table  1  we  list  simulation  estimates,  numerical  values  from  the  exact  formula  and  from 
the  light-heavy  traffic  interpolation.  A  comparison  of  the  exact  values  to  the  light-heavy 
traffic  interpolation  shows  that,  in  this  case,  the  agreement  is  excellent.  Since  we  expect 
the  approximation  to  be  asymptotically  exact  at  the  endpoints  p  =  0  and  p  —  1,  it  is  not 
surprising  that  the  largest  errors  occur  in  moderate  traffic. 


Tail  probability  P  p  [q  >  4] 

P 

Simulation 

Approximation 

Error  (%) 

0.1 

1.1271e-04±1.7% 

9.0718e-05 

19.51 

0.2 

1.3444e-03±1.7% 

9.4753e-04 

29.52 

0.3 

6.1736e-03±0.9% 

5.9952e-03 

2.89 

0.4 

1.9246e-02±0.6% 

1.4415e-02 

25.10 

0.5 

4.7745e-02±0.5% 

3.9894e-02 

16.44 

0.6 

1.0292e-01±0.4% 

9.5777e-02 

6.94 

0.7 

2.0035e-01±0.3% 

1.9757e-01 

1.39 

0.8 

3.6093e-01±0.3% 

3.6379e-01 

-0.79 

0.9 

6.1407e-01±0.3% 

6.1902e-01 

-0.81 

Table  2:  P  p  [q  >  4]  for  deterministic  burst  duration  cr  =  3. 

In  the  same  setup,  we  next  consider  the  tail  probability  P  f)  [q  >  4].  From  Table  2  we 
see  that  although  the  approximation  yields  estimates  in  the  correct  order  of  magnitude, 
the  errors  are  substantial  when  not  in  the  moderate-to-heavy  traffic  regime.  This  can  be 
explained  as  follows:  When  a  =  3,  in  order  for  the  queue  to  build  up  to  4  at  least  3  sources 
should  be  simultaneously  active.  Note  that  the  light  traffic  component  of  the  approximation 
consists  of  the  second  derivative,  which  can  be  obtained  by  considering  sample  paths  with 
at  most  two  source  arrivals  in  the  system.  Thus,  any  effects  due  to  the  activation  of  more 
than  two  sources  are  not  adequately  accounted  for  in  light  traffic. 
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Uniform  Specializing  (18)  to  the  case  where  a  is  uniformly  distributed,  P  [a  —  n]  =  1/M, 
n  —  1,  2, . . . ,  M,  yields 


P  P  [q  >  b\ 


1  [M  >  (i  -  p)(,]  (l  +  5(M 

x  (M  -  (1  -  p)b )2  +  p3  exp  f-6 


(i-pW2) 

6  =  0,1,.... 


Tail  probability  Pp  [q  >  b] 

Buffer  size 

Simulation 

Approximation 

Error  (%) 

0 

4.4907e-02±0.2% 

4.5333e-02 

-0.95 

2 

1.1747e-02±0.5% 

1.1399e-02 

2.97 

5 

1.4543e-03±1.4% 

9.7380e-04 

33.04 

8 

1.9456e-04±3.7% 

2.4379e-04 

-25.30 

10 

5.1620e-05±7.0% 

1.0186e-04 

-97.31 

Table  3:  Traffic  intensity  p  =  0.2;  a  ~  nniform(l,  5). 

For  M  =  5  we  compare  simulation  vs  approximation  in  Tables  3  and  4,  for  traffic  in¬ 
tensities  p  =  0.2  and  p  =  0.8  respectively.  Once  more,  the  approximation  is  very  sharp 
for  small  buffer  sizes.  As  the  buffer  size  increases  beyond  the  maximum  burst  length  and 
the  true  probabilities  become  smaller,  the  approximation  lingers  on  in  the  correct  order  of 
magnitude,  but  it  clearly  deteriorates  away  from  heavy  traffic.  Eventually,  as  the  buffer  size 
tends  to  infinity,  the  interpolation  approximation  overestimates  the  actual  probabilities. 


Tail  probability  Pp  [q  >  b] 

Buffer  size 

Simulation 

Approximation 

Error  (%) 

0 

6.5489e-01±0.1% 

6.6133e-01 

-0.98 

10 

2.0086e-01±0.4% 

1.9161e-01 

4.60 

20 

6.3303e-02±0.9% 

5.8057e-02 

8.29 

30 

1.9964e-02±1.7% 

1.9406e-02 

2.79 

40 

6.2827e-03±3.1% 

6.5188e-03 

-3.75 

50 

1.9703e-03±5.6% 

2.1897e-03 

-11.13 

Table  4:  Traffic  intensity  p  =  0.8;  a  ~  uniform(l,  5). 


Geometric  Taking  a  to  follow  a  geometric  distribution,  P  [cr  >  n]  =  pn,  n  =  0, 1, . . .,  we 
obtain  from  (18)  that 

ppfe  >  b)  «  y(l  -p)(l  +pfp2{l~p)h  +  p3exp  ^-2  Y“^(1  ~  P)b)  ’  b  =  0,l,... 
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Figure  1:  Geometric  p  =  0.8  burst  duration. 


As  an  example  we  set  p  =  0.8  and  plot  simulated  and  approximate  values  in  Figure  1,  for 
traffic  intensities  p  =  0.1,  0.4,  0.6  and  0.9.  In  all  cases  confidence  interval  widths  were  within 
10%  of  the  mean  and  are  not  shown.  The  linear  decrease  of  the  simulated  values  suggests 
an  exponential  decay  of  the  queue  length  distribution,  in  agreement  with  large  deviations 
results.  Figure  1  clearly  indicates  that  the  heavy-light  traffic  interpolation  is  sufficient  for 
providing  rough  estimates  for  a  wide  range  of  probabilities  and  buffer  sizes. 


6  Extensions 


It  is  apparent  from  the  developments  of  Section  3  that  the  light  traffic  results,  as  stated  in 
Proposition  (3.1),  do  not  cover  several  interesting  distributions  belonging  to  the  subexpo¬ 
nential  family.  Such  is  for  example  the  lognormal  distribution,  which  violates  Assumption 
(A)  despite  having  finite  kth  moment  for  every  k  =  0, 1, . . ..  It  is  natural  to  expect  that 


a 


Assumption  (A)  can  be  relaxed  to  require  that  E 
order  for  Proposition  (3.1)  to  go  through.  This  won 
dependence,  for  which  E  [a2]  =  oo.  We  present  some  heuristics  below. 


be  finite,  for  appropriate  k  >  2,  in 
d  still  not  address  the  case  of  long-range 


Long-range  dependence  For  the  Pareto  distribution  P  [a  >  n]  =  n_",  n  =  1,2,..., 
1  <  a  <  2,  not  only  Assumption  (A)  fails,  but,  as  mentioned  in  Section  3,  (6)  yields 
infinity.  This  suggests  that  Pf,  [q  >  6]  may  not  be  an  analytic  function  of  p  under  long-range 
dependence.  When  c  =  1  simulation  work  in  light  traffic  does  not  preclude  the  possibility  that 
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limp  "P p  [q  >  b]  is  the  sought  after  lion-trivial  limit.  On  the  other  hand,  the  heavy  traffic 

result  of  Proposition  3.2(b)  hints  at  developing  an  approximation  around  the  normalized  rv 
(1  —  p)1/!"-1)  q.  These  considerations  lead  us  to  postulate  that,  when  c  =  1, 


limp  "P,  [q  >  b]  -  I\(b), 


(21) 


for  some  unknown  mapping  K  :  R+  — >  R+.  Then,  taking  advantage  of  Proposition  3.2(b) 
we  propose  the  approximant 


Pp  [q  >  b] 


(a  -  1)E  [a]  1  -p  j\ 
P(2  -  a)  p«  )' 


c  =  1. 


(22) 


This  expression  is  in  agreement  with  the  heavy  traffic  limit  (8).  In  addition,  from  the 
Mi t tag  -LefTler  function  asymptotics  given  in  [5,  p.  207],  we  have 


Ea-i(-x) 


1  1 

x  T(2  —  a ) 


(x  oo ) 


which  ensures  that,  as  p  — »  0,  approximation  (22)  conforms  with  the  conjectured  fight  traffic 
limit  (21). 


Assessing  the  performance  of  (22)  requires  numerical  evaluation  of  the  Mittag-Leffler 
function.  In  general,  calculation  based  on  the  series  expansion  (9)  is  not  recommended. 
Instead,  one  can  invert  the  Laplace  transform  of  the  Mittag-Leffler  law  by  contour  integration 
along  a  suitably  chosen  path  in  the  complex  plane.  Doing  so  [15]  we  arrive  at  the  alternate 
expression 


sin(z/7r)  rn/2  e  0rtan6')1/l'  ^ 

isit  Jo  1  +  sin(20)  cos(z/7r) 


x  >  0,  0  <  v  <  1, 


which  is  evaluated  by  numerical  integration.  We  then  test  approximation  (22)  for  traffic 
intensities  p  =  0.2,  0.5  and  0.8.  Under  long-range  dependence  simulation  estimates  converge 
very  slowly;  moreover  confidence  intervals  based  on  the  regenerative  method  cannot  be 
constructed,  because  the  underlying  period  has  infinite  variance.  In  the  results  shown  the 
runs  were  109  time  slots  long,  and  by  that  time  the  estimates  had  stabilized.  The  log-log 
scale  plots  in  Figures  2  and  3  correspond  to  two  Pareto  distributions  with  parameters  a  =  1.5 
and  a  =  1.7.  Observe  that  the  heavier  a  =  1.5  Pareto  tail  induces  larger  tail  probabilities 
than  a  =  1.7,  at  the  same  traffic  intensities.  In  both  Figures  2  and  3  we  see  that  simulated 
and  approximate  values  are  very  close,  suggesting  that  expression  (22)  provides  a  satisfactory 
approximation.  Note  also  the  almost  linear  shape  of  the  curves  in  log-log  scale,  reflecting 
the  power  law  asymptotics  of  the  queue  size  distribution  announced  in  [7,  8,  9]. 
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